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A novel investigation of the nature of intermittency in incompressible, homogeneous and isotropic 
turbulence is performed by a numerical study of the Navier-Stokes equations constrained on a 
fractal Fourier set. The robustness of the energy transfer and of the vortex stretching mechanisms is 
tested by changing the fractal dimension, D, from the original three dimensional case to a strongly 
decimated system with D = 2.5, where only about 3% of the Fourier modes interact. This is a 
unique methodology to probe the statistical properties of the turbulent energy cascade, without 
breaking any of the original symmetries of the equations. While the direct energy cascade persists, 
deviations from the Kolmogorov scaling are observed in the kinetic energy spectra. A model in terms 
of a correction with a linear dependency on the codimension of the fractal set, E{k) ~ , 

explains the results. At small scales, the intermittency of the vorticity field is observed to be 
quasisingular as a function of the fractal mode reduction, leading to an almost Gaussian statistics 
already at D ~ 2.98. These effects must be connected to a genuine modification in the triad-to-triad 
nonlinear energy transfer mechanism. 

PACS numbers: 


Understanding and controlling the energy transfer pro¬ 
cess in turbulent flows is a key problem for a broad 
range of fields, such as astrophysics [T], atmospheric 
or ocean sciences [5], mathematics and engineering [3]. 
The main obstacle hampering theoretical, numerical, and 
phenomenological advancements is intermittency: the 
presence of strong non-Guassian and out-of-equilibrium 
velocity fluctuations in a wide range of scales [IHH]- The 
energy cascade is fully chaotic, nonlinear and driven by 
the vortex-stretching mechanism, i.e., the tendency of 
the flow to amplify vorticity in thin, long filaments. A 
long debate exists whether or not the presence of such 
geometrical structures is correlated to the non-Gaussian 
statistics observed at small scales [3]. Many authors 
have focused on a vortex-by-vortex analysis, looking for 
the signatures of quasisingularities or extreme events 
due to specific dynamical properties of the Navier-Stokes 
(NS) equations [9HT5]. Other approaches are based 
on a Fourier description, such as closure [M HU and 
renormalization-group theories [TSH^O] . 

In this Letter, we investigate the origin of intermittency 
in turbulent flows by a novel strategy. The idea consists 
of modifying the nonlinear interactions of Navier-Stokes 
equations, without introducing extra forces and without 
breaking the symmetries of the original equations, such 
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as statistical homogeneity, isotropy and rescaling prop¬ 
erties in the inviscid limit. To achieve this, we adopt 
a recently proposed numerical methodology m to solve 
NS equations on a preselected, multiscale (fractal) set 
of Fourier modes. This allows us to control the number 
of degrees of freedom participating in the non-linear dy¬ 
namics by changing one free parameter only, the fractal 
dimension of the Fourier set, D. For D = 3, the origi¬ 
nal problem is recovered. In the sequel, we describe the 
methodology, the numerical setup and the main results 
concerning the quasisingular effect of fractal mode reduc¬ 
tion on the small-scale intermittency. 

Methodology. Fractal mode reduction is realized via the 
Fourier decimation operator V^, acting in the space of 
divergece-free velocity fields as follows m- We define 
v(x, t) and u(k, <) as the real and Fourier space repre¬ 
sentation of the velocity field m D — 3, respectively. The 
decimated field, v^(x, t), is obtained as: 

v^(x,t) ='P^v(x,t) = ^ 7 ku(k, t). (1) 

The random numbers 7 k are quenched in time and are: 

1, with probability hk , 

0, with probability 1 — hk, A: = |k|. 

The choice for the probability hk oc {k/kQ)^~^, with 
0 < D < 3 ensures that the dynamics is isotropically 
decimated to a Z?-dimensional Fourier space. The fac- 
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tors hk are chosen independently and preserve Hermi- 
tian symmetry 7 ^ = 7 _fe so that is self-adjoint. The 
NS equations for the velocity field decimated on a fractal 
Fourier set are then defined as: 

dtv^ + + . (3) 

At each iteration the nonlinear term, A(v,v) = —v • 
Vv — Vp, is projected on the quenched fractal set, to 
constrain the dynamical evolution to evolve on the same 
Fourier skeleton at all times. Similarly, the initial condi¬ 
tion and the external forcing must have a support on the 
same set of Fourier modes. Let us notice that the pro¬ 
jection acts as a self-similar Galerkin truncation. In the 
(L^) norm, ||v|| cx / |v(x)pd^a:, the self-adjoint operator 
commutes with the gradient and viscous operators. 
Since = v^, it then follows that both energy and 

helicity are conserved in the inviscid and unforced limit, 
exactly as in the original problem with D = 3. As a result 
of the Fourier decimation, the velocity field is embedded 
in a three dimensional space, but effectively possesses a 
number of Fourier modes that grows slower with decreas¬ 
ing D. The degrees of freedom inside a sphere of radius 
k go as ^dof{k) ~ k°. 

This idea, introduced for the first time in m, has 
been used to test the hypothesis that two-dimensional 
turbulence in the inverse energy cascade approaches a 
quasiequilibrium state [52]. Fourier decimation methods 
are not new for hydrodynamics: we mention protocols 
with a specihc degree of mode reduction [55H55] , and the 
extreme truncation criterion of shell models for the tur¬ 
bulent energy cascade |26|. Results are puzzling. For NS 
equations at small Reynolds numbers |23j . intermittency 
strongly depends on the amount of scales resolved in the 
inertial range. However, in the case of shell models, in¬ 
termittency is observed to be a function of the effective 
dimension m, and coinciding with the one measured 
in the original NS equations |28j . when energy and he¬ 
licity are the two inviscid invariants. Note that fractal 
Fourier mode-reduction changes also the relative popu¬ 
lation of local-to-nonlocal triadic interactions |S], since 
triads with all modes in the high wave number range 
have a larger probability to be decimated. Furthermore, 
being an exquisitely dynamical approach, it is different 
from a posteriori filtering techniques, largely exploited to 
analyze turbulent data [55] . 

A pseudospectral method is adopted to solve Eqs. ^ 
in a periodic box of length L = 27r at resolution N = 
1024^ and 2048^, dealiased with the two-thirds rule; time 
stepping is implemented with a second-order Adams- 
Bashforth scheme. A large-scale forcing [5D] keeps the 
total kinetic energy constant in a range of shells, 0.7 < 
|k| < 1.7, leading to a steady, homogeneous and isotropic 
flow. We performed direct numerical simulation (DNS) 
runs by changing the fractal dimension 2.5 < D < 3, the 
spatial resolution, the viscosity and the realization of the 
fractal, quenched mask. Table [^ summarizes the relevant 


D 

3 

2.999 

2.99 

2.99 

2.98 

2.98 

2.8 

2.5 

N 

1024 

1024 

1024 2048 

1024 2048 

1024 

1024 

Mr 

100% 99% 

93% 

92% 

87% 

85% 

25% 

3% 

V 

0.75 

0.75 

0.95 

0.70 

0.75 

0.70 

0.90 

0.65 

Mt 

10 

10 

11 

10 

11 

7 

10 

20 

E 

3.1 

3.1 

3.1 

3.3 

3.3 

3.5 

4.1 

5.4 

V- lO'* 

6.0 

6.0 

6.0 

2.0 

6.0 

2.0 

6.0 

1.5 

Re ■ 10"® 

3.9 

3.9 

3.8 

11.8 

3.9 

12.1 

4.0 

15.4 


TABLE I: DNS parameters. D-. fractal dimension; A: grid 
resolution; Mr', percentage of surviving Fourier modes; 7 : 
Kolmogorov length scale in unit of the grid spacing Ax = 
L/N', Nt'. number of large-scale eddy-turnover-times in the 
steady state; E-. total kinetic energy (u^)/2; v. viscosity; Re: 
Reynolds number Re = E^^"^Lfv, where the integral scale L 
is estimated from the kinetic energy spectrum. 


parameters. 

Results. The starting point of our analysis is the shell- 
to-shell energy transfer in the Fourier space. Following 
the notation adopted in Ref. [S], we write the energy 
spectrum for a generic flow in dimension D as: 

E^{k)= [ d^fci7ki /'d^A:27k2(u(ki)u(k2)) , (4) 
J |ki \—k J 

where the decimation factor 7 k takes into account that 
the Fourier mode k is active with probability hk. Sim¬ 
ilarly, we can write for the energy flux across a Fourier 
mode fc, n^(/c) = J^^^^^^i^d^kidtE{ki): 

n^(fc) = [ d^fci7ki f d^fc2d^A:37k2 7k35'(ki|k2,k3), 
-'|ki|<fc J 

_ ( 5 ) 

where the explicit form of the symmetric triadic cor¬ 
relation function is [5T]: S'(ki|k 2 ,k 3 ) = — /TO[((ki • 

u(k 3 ))(u(ki) • u(k 2 ))) -t- ((ki • u(k 2 ))(u(ki) • u(k 3 )))]. 
Supposing a power-law behavior of the velocity fluctua¬ 
tions u{k) ~ k~°‘, we can predict a self-similar scaling 
of the energy flux as n'°(Afc) ~ In this 

expression, the rescaling factor A^^ is due to the integral 
over the variables (ki,k 2 ,k 3 ), while A^“^“ comes from 
the triadic nonlinear term. 

If a constant energy flux develops in the inertial range of 
scales, the following dimensional relation holds: 

a = D+ 1/3E^{k)r~.^k^-^ E^^Hk), ( 6 ) 

where E^'^^fk) ~ k~^/^ is the Kolmogorov 1941 spec¬ 
trum expected for the original case with D = 3. In 
the previous dimensional argument, the tiny inter¬ 
mittent corrections to the spectrum scaling exponent 
are neglected [35], while prefactors are omitted for 
simplicity. The relation (|^ is obtained by notic¬ 
ing that, because of homogeneity, we have that 
( u(ki) • u(k 2 )) oc F(ki)5^(ki -f k 2 ), and by also noticing 
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FIG. 1: (Top): Log-Log plot of the mean kinetic energy spec¬ 
tra for different D; in the inset, the mean kinetic energy fluxes. 
(Bottom): Compensated energy spectra (k) -3+^ vs 

the wave number k; in the inset, the compensation is done 
with the Kolmogorov prediction, (k) k^^^. 


that the decimation projector verifies the identity 
(7k)^ = 7k- As a result, the dynamical effect of Fourier 
decimation is to make the energy spectrum shallower 
than the Kolmogorov one for three-dimensional turbu¬ 
lence, also predicting the existence of a critical dimension 
D = 7/3, when the spectrum becomes ultraviolet diver¬ 
gent m in the limit of zero viscosity. By decreasing D 
in the presence of a forward energy cascade, the system 
has fewer modes available to transfer the same amount 
of energy (see Table [l| , and the velocity field becomes 
increasingly rougher. 

In the upper panel of Fig. we plot the energy spectra 
and the associated energy fluxes, for all DNS. It shows 
that when increasing the grid resolution for fixed D, 
from N = 1024 to TV = 2048, no appreciable differences 
are observed, indicating that the presence of a forward 
energy cascade appears robust and Reynolds indepen¬ 
dent. In the lower panel we also show that the spectra 
compensate well with the prediction ([^, while they 
fail to satisfactorily compensate with the classical K41 
prediction when D < 3. We stress that the latter result 
is significant for D = 2.5 and D = 2.8 only. For the 
other dimensions D, the effect is so small that it might 



FIG. 2: Probability density function of the vorticity compo¬ 
nent uix, normalized to its standard deviation. Data refer to 
simulations at resolution N = 1024. In the inset, mean square 
vorticity {uP') versus the fractal dimension deficit, 3 — D. 


fall within the intermittent correction of the original 
Navier-Stokes case at D = 3. Moreover, the effect 
of the quenched disorder is robust: spectra obtained 
with different realizations of the mask do not show any 
statistically significant difference (not shown). 

Figure]^ (upper inset) shows that, when decreasing the 
fractal dimension D, the mean energy transfer towards 
small scales is almost unchanged, i.e. the hypothesis 
leading to the relation (§ is well verified. On the other 
hand, temporal fluctuations of the energy flux increase 
with D (not shown). 

It might be argued that the effect of fractal Fourier 
decimation is purely geometrical and that the main 
dynamical processes are unchanged. To show that this is 
not the case, it is useful to analyze the effect of a static 
Fourier decimation. This can be done by considering 
snapshots of D = 3 turbulence, and applying the fractal 
decimation as an a posteriori filter. It is immediate to 
realize that the effect of the static decimation on the 
spectrum is E^{k) ~ implying that the 

geometrical action of the decimation goes in the opposite 
direction of the dynamical one. 

We now consider the dynamical effect of the fractal 
Fourier decimation on the small-scale structures, by 
focusing on the statistics of the vorticity field in the 
real space. In Fig. we plot the probability density 
function (PDF) of the vorticity field, normalized with its 
standard deviation. We note that already at D = 2.99, 
vorticity fluctuations have changed their intensity of 
I order of magnitude, despite the fact that the mean 
enstrophy is practically unchanged. Even more re¬ 
markable, intermittent fluctuations disappear already 
at D = 2.8, where a quasi-Gaussian vorticity PDF 
is measured. The transition towards Gaussianity is 
better quantified considering the vorticity kurtosis. In 
Fig.d we compare results of the fractally decimated 
NS equations, with those obtained from the application 
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3-D 

FIG. 3: Lin-log plot of vorticity kurtosis vs the dimension 
deficit 3 — D. The upper curve is obtained from the applica¬ 
tion of the a posteriori static mask on D = 3 velocity field 
snapshots. The lower curve comes from the Fourier decimated 
DNS. Data for D = 2.99 and D = 2.5, with a different real¬ 
ization of the fractal projection, are indistinguishable. Data 
at resolution N = 2048 are also plotted. 


of the a posteriori static mask on three-dimensional 
turbulence. The dynamical decimation makes a very fast 
transition towards a Gaussian behaviour, such that at 
D = 2.98 the kurtosis has decreased by 30%, to already 
approach the Gaussian value at Z? = 2.8. In the case 
of thea posteriori static decimation, vorticity kurtosis 
assumes the Gaussian value only at 15 = 2.5, while 
staying almost unchanged in the range D > 2.98. Such 
a strong difference clearly indicates that constraining 
the dynamics to a subset of modes is critical for the 
complete development of intermittency in real space. 

In Fig.|^ we plot the kurtosis of the longitudinal velocity 
increment, K{6vr) = ((^rV • r)"‘)/((5rV • r)^)^. Notice the 
sharp transition leading to an almost scale-independent, 
nonintermittent behavior for Z5 < 3. Moreover, at fixed 
fractal dimension, the effects of increasing the Reynolds 
number is to further reduce the intermittent corrections. 
Conclusions. Fractal mode-reduction is a new route 
to design numerical simulations to tackle the problem 
of intermittency and to potentially develop multiscale 
models of turbulence. The first non trivial result is the 
robustness of the energy flux under modereduction. An 
inertial range of scales with a constant-flux solution is 
observed when D is changed and few Fourier modes sur¬ 
vive, at least in the parameter range investigated here. 
This is in agreement with the observation that Galerkin 
truncations do not alter the inviscid conservation of 
quadratic quantities, preserving the existence of exact 
scaling solutions for suitable third-order correlation 
functions (see appendix of Ref. [34]). The Fourier 
spectrum gets a power law correction that can be 
predicted by a dimensional argument. Second and more 
striking, small-scale intermittency is quickly reduced for 
D < 3 and it almost vanishes already at 15 = 2.98. As 



FIG. 4: Lin-log plot of longitudinal velocity increment kur¬ 
tosis K{Svr) vs the scale separation normalised by the Kol¬ 
mogorov scale, r/rj. For all data sets, error bars are within 
the symbol width. 


a consequence, the presence or absence of some of the 
Fourier modes strongly modify the fluctuations of all the 
others, suggesting the possibility that intermittency is 
the result of percolating dynamical properties across the 
whole Fourier lattice |55] . 

Because of the spectrum modification, the scaling expo¬ 
nent of the second order longitudinal structure function 
becomes (2 + {D — 3), where (2 is the measured in the 
D — 3 homogeneous and isotropic case. This observation 
would suggest that, for the dimension deficit 3 — Z5 < 1, 
one may obtain corrections to the scaling exponents 
proportional to 3 — D, and the anomalous exponents 
might be computed perturbatively in the dimension 
deficit. If this is the case, the critical dimension D,, 
is estimated as the value of D where the dimensional 
scaling is recovered, namely C 2 + Z5c — 3 = 2/3. It 
gives Dc ^ 2.96 not far from the value of D at which 
intermittency is observed to vanish in the DNS. However, 
there is no reason to assume that anomalous exponents 
can be computed perturbatively in 3 — 15. In fact, as 
mentioned above, intermittency might well be the result 
of multiscale interactions in Fourier space, needing all 
degrees of freedom to develop. Hence, in the presence of 
even a tiny decimation, NS singular solutions responsible 
for the anomalous scaling no longer exist. This would 
also explain why we observe. Fig. [^ a reduction of 
intermittency by increasing the resolution at fixed 
D. Additionally, phenomenological cascade models |4] 
would be unable to explain the results, as well. 

In the light of our results, Fourier decimation can also 
be seen as a way to introduce a control parameter and 
modify the scaling properties of the system, similarly 
to what happens for NS equations stirred by a random, 
power-law forcing [35H^ . In the latter case, perturba¬ 
tive or semianalytic calculations gniiiT] give indications 
on the reasons why anomalous corrections should cancel 
out for specific values of the control pararameter. Also, 
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in Ref. |37j it is numerically shown that when the random 
injection becomes the dominant scaling contribution in 
the inertial range, a transition to a Gaussian statistics 
is observed for the velocity increments. In the present 
case, however, the connections between the observed 
transition to a Gaussian behavior, and a possible 
renormalised perturbation theory are to be explored. 
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